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Improving the accuracy of forecast models for physical systems such as the atmosphere is a 
crucial ongoing effort. Errors in state estimation for these often highly nonlinear systems has been 
the primary focus of recent research, but as that error has been successfully diminished, the role 
of model error in forecast uncertainty has duly increased. The present study is an investigation 
of a particular empirical correction procedure that is of special interest because it considers the 
model a "black box", and therefore can be applied widely with little modification. The procedure 
involves the comparison of short model forecasts with a reference "truth" system during a training 
period in order to calculate systematic (1) state-independent model bias and (2) state-dependent 
error patterns. An estimate of the likelihood of the latter error component is computed from the 
current state at every timestep of model integration. The effectiveness of this technique is explored 
in two experiments: (1) a perfect model scenario, in which models have the same structure and 
dynamics as the true system, differing only in parameter values; and (2) a more realistic scenario, 
in which models are structurally different (in dynamics, dimension, and parameterization) from the 
target system. In each case, the results suggest that the correction procedure is more effective for 
reducing error and prolonging forecast usefulness than parameter tuning. However, the cost of this 
increase in average forecast accuracy is the creation of substantial qualitative differences between 
the dynamics of the corrected model and the true system. A method to mitigate the structural 
damage caused by empirical correction and further increase forecast accuracy is presented. 
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INTRODUCTION 



Advances in computational power and increasingly ac- 
curate techniques for estimating the current state of the 
Earth's atmosphere have significantly improved numeri- 
cal weather prediction (NWP) [T]-[3]. As state estimation 
error is reduced due to improved methods of data assim- 
ilation, error in the model tendency plays an increasing 
role in the uncertainty of predictions at every temporal 
and physical scale [4]-[6]. 

In 1978, Leith introduced a statistical technique to cor- 
rect model tendency error, in which short model forecasts 
are compared to a time series of reference "truth" to 
estimate both state-independent model bias, and state- 
dependent error components which are approximated by 
a least-squares linear function of the model state [?]• 
More recently, empirical correction has been employed 
with success in atmospheric models with relatively few 
degrees of freedom (e.g. N = O(IO^) in [8 ), and low- 
dimensional modifications of the technique involving, 
for example, singular value decomposition (SVD) of the 
state-dependent correction operator (the least-squares 
linear function proposed by Leith) have proven successful 
in models with as many has O(IO^) degrees of freedom 
[9]. In this study, we apply the original technique devel- 
oped by Leith to simple three-dimensional Lorenz models 
[To] , where in addition to testing the effectiveness of em- 
pirical correction, we aim to understand the dynamical 
ramifications of a statistical approach to the correction 
of model tendency error. 

The model tendency M(x) is defined as the change in 
state variables over one timestep of numerical integra- 



tion, which we denote as a time derivative: 

M(x) = x^ (1) 

where x is the atmospheric state-vector, typically with 
0(10^^) degrees of freedom for NWP. Note that x^ rep- 
resents the model state, whereas x-^ will represent the 
true system state, in terms of model variables. 

Given the state x of a physical system, M consists of 
all the known physics, forcings, and parameterizations of 
sub grid-scale processes. To make a one-timestep model 
forecast, we approximate the true change in state vari- 
ables over that time by the model tendency, x-^ ~ ^(x), 
and the error 6M in that approximation is called the 
tendency error: 



SM 



M(x) 



(2) 



Even with perfect estimates of the current state of the 
atmosphere, the model tendency error would quickly sep- 
arate forecasts from the truth, due to the atmosphere's 
chaotic dynamics [11 . As a result, techniques for reduc- 
ing the model tendency error represent a current and ac- 
tive research area, and those that are applicable indepen- 
dent of the specific model are of special interest. Clearly, 
one would like to improve the physics represented by M 
from first principles. In what follows, we assume this im- 
provement has met the limit of diminishing returns, and 
move towards a statistical approach. 

The general strategy of empirical correction is to com- 
pare short forecasts generated by a model to observations 
of the system being modeled over some training period. 
If the state-space of the system is well represented by the 
training period, and the model is a reasonable approx- 
imation of the true system, the forecast error statistics 



can be used to create an empirical correction that pushes 
the model closer to the truth at each timestep of numer- 
ical integration. Adjusting the model every timestep re- 
duces the nonlinear growth of tendency error, providing 
more effective error reduction than a posteriori statis- 
tical correction [12 . This strategy is similar to nudg- 
ing or Newtonian relaxation in a data assimilation (DA) 
context, where one is assimilating observations, except 
that here we are nudging with predicted, rather than ob- 
served, forecast error. 

The present study is an investigation of a three-step 
empirical correction procedure inspired by the work of 
Leith [7 , DelSole and Hou [8 , and more recently Dan- 
forth et. al. 0[T3]. We first test its effectiveness in syn- 
chronizing Lorenz systems [10 with varied parameter- 
values in a perfect model scenario. We then apply the 
correction to an alternative model derived by Ehrhard 
and Miiller [14 tuned to approximate the evolution of 
a toroidal thermosyphon, an experimental analogue to 
the original Lorenz system. The "true" climate is repre- 
sented by a long-time, high-dimensional computational 
fluid dynamics (CFD) simulation of the thermosyphon. 
An analysis^ which is an approximation of the true sys- 
tem state in terms of model variables, is then created by 
three-dimensional variational (3DVar) data assimilation 
and used for training and verification of the empirical cor- 
rection. This process mimics the application of empirical 
correction in an operational NWP setting. We also verify 
the corrected model by direct comparison with observa- 
tions of the truth. 

The results in each experiment suggest that the cor- 
rection procedure is effective for reducing error. How- 
ever, there is an associated cost of this short-term error 
reduction, which is evidenced by substantial qualitative 
differences between the dynamics of the corrected model 
and the true system, differences that were not present in 
the uncorrected model. Introduction of system-specific 
knowledge into the correction procedure is shown to mit- 
igate some of that cost, while also improving error statis- 
tics further than the entirely general procedure. 

The paper is structured as follows. In Sec.|n]we define 
the procedure for a general model M. The application of 
the technique in the perfect model scenario is addressed 
in Sec. |III[ and in Sec. |IV| we present the thermosyphon 
model correction. Finally, we discuss the results and con- 
clude the paper in Sec. [Vl 



II. EMPIRICAL CORRECTION 

The correction procedure employed in this experi- 
ment consists of three steps: (A) training, (B) state- 
independent correction and (C) state-dependent correc- 
tion. The state-independent correction can be thought of 
as aligning the time-average of the model state with that 
of the true state. Likewise, the state-dependent correc- 
tion can be considered an alignment of the model vari- 
ance with that of the truth. To determine the correction 



terms, we compare short model forecasts to observations 
of the true system over a training period in a process 
called direct insertion. 



A. Training 

In general, comparing model forecasts to a true phys- 
ical system requires estimates of the true system state 
in terms of the model state- variables. Consider a vector 
time-series x-^(t) of such estimates, which we will call the 
"reference truth" . The amount of time /i, measured in 
model timesteps, between estimates is called the analysis 
window; we assume it to be constant. 

The process of direct insertion begins with the gen- 
eration of a time-series x^(t) of duration-/^ model fore- 
casts, where each forecast in the time-series is initialized 
from the previous state in the reference truth. The first 
vector in the series, for example, will be x^(to + h)^ 
which is the model state resulting from an /i-timestep 
forecast started with initial condition x^(to). Subtract- 
ing each of the model forecast states x^ (t) from the cor- 
responding reference true state x-^(t) produces a third 
time-series Ax(t) which represents the forecast errors af- 
ter h timesteps. These errors result from differences be- 
tween the model rate of change for each variable and the 
true rate of change, and they are commonly referred to 
as analysis corrections (or increments). See Fig. [l] for a 
schematic of the procedure. 
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FIG. 1. The direct insertion procedure for comparing short 
model forecasts to the truth to obtain a time series of analysis 
corrections, x^ represents a time series of the reference truth, 
and the analysis window represents the number of timesteps 
between estimations of the true system state, x^ represents 
a time series of forecasts with duration equal to the analy- 
sis window, each of which is started from the previous true 
state. The time- average of the analysis corrections (Ax) di- 
vided by the number of timesteps in the analysis window h 
approximates the average (state-independent) model bias b. 



Finally, we separate each of the time series into anoma- 



lous and time-average components: 



x^'(t) 
Ax'(i) 



= x^(t) - (x^) 



M 



(t)-{x-) 
Ax(t) - (Ax) 



(3) 



where the expectation operator (•) denotes averaging over 
the training period, and the primes denote anomalies^ 
which are differences from the mean. The time- average 
components wih be used for state-independent correction 
as described in Sec. |IIB| and the anomaly time-series wih 
be used for state-dependent correction as detailed in Sec. 

inn 



B. State-independent correction 

We turn our attention first to a state-independent cor- 
rection of the form 



M*(x) = M(x)+b 



(4) 



where the constant vector b is the average model error 
(bias) to be determined. 

Recall that our goal here is to empirically align the 
time-averages of the state- variables in the model to those 
of the true system. We call the time-averaged true system 
state the climatology^ and we approximate it by (x-^), 
the average of the reference true state over the entire 
training period. The average of the analysis corrections 
(Ax) over the training period provides an estimate for 
the systematic, state-independent error generated by the 
model during the analysis window, as explained in Fig. 
[l] Dividing by the number of timesteps in the analysis 
window, then, we approximate the model bias by b = 
(Ax)//i, and the bias-corrected model tendency is thus 
given by 



M*(x) = M(x) 



(Ax) 
h 



(5) 



Note that at this point we are approximating the model 
tendency error SM by the model bias b alone. We also 
wish to estimate any component of error that may depend 
on the system state, by approximating SM ^ h -\- Lx\ 
where L is a matrix operator to be described in the next 
section. 



C. State-dependent correction 

To generate a linear state-dependent correction oper- 
ator L, we follow Leith [7 and DelSole and Hou [8 , by 
first recomputing the forecast and correction time series 
in Fig. [l] using the state-independent corrected model, 
M*, and then decomposing into mean and anomalous 
components. We seek an improved model of the form 



where M+ includes both stages of correction. Letting 
g = 6M^ = x^ — [M* (x) + Lx'] be the tendency error 
of the improved model, we minimize the expected square 
tendency error (g^^g), (where g^^ is the transpose of g), 
with respect to L. The minimization results in the for- 
mula 



L = (Ax'x^^^)(x 



.r^rtry 



1 _ 



^Axx 



tC: 



-1 



(7) 



as explained by Danforth et. al. [9 . C^xx^ is the aver- 
age over the training sample of cross covariance matrices 
obtained by taking the outer product [Ax^(t)] • [x^ {t)Y^ 
for each time t. C^^t^^t is the average true-state covari- 
ance matrix, and L is known as Leith's state-dependent 
correction operator. 

When L operates on the current anomalous state x', 
we can think of it as doing two things: (1) C~^ ^ op- 
erates on x', effectively relating the current state to the 
reference truth in the dependent sample, i.e. giving the 
best representation of the current state in terms of past 
states; and then (2) C^^^x^ operates on the result, deter- 
mining what correction should be made. This allows the 
model correction to adjust to different regions of state- 
space, and explains why the state-dependent correction 
can be thought of as attempting to align the model and 
true-state variances. 

In the next section, we describe the application of this 
three-stage procedure to align Lorenz systems with dif- 
ferent parameter values in a perfect model scenario, and 



in Sec. IV we describe its application to couple a low di- 
mensional model to a high dimensional toy climate simu- 
lation. We also note here that the term corrected will im- 
ply the application of both state-independent and state- 
dependent correction, unless explicitly stated otherwise. 



III. PERFECT MODEL SCENARIO 

As a first step in the investigation of the correction 
technique, we consider its application to a model origi- 
nally studied by Lorenz [10]. The system of equations 
Q represents fluid flow between two plates, Rayleigh- 
Benard convection, in which convection cells form for 
certain parameter ranges. However, with only slight 
modification (the details of which appear in Sec. IV A), 
they also describe the flow in a natural convection loop 
[Ml [15] . Lorenz systems are covered exhaustively in pub- 
lication [16l dl] , and thus provide a familiar platform on 
which to perform preliminary tests of strategies for pre- 
dicting the future state of chaotic systems. 



x^ 



M+(x) = M*(x) + Lx' 



(6) 



A. Experimental design 

In this perfect model scenario, the true system and 
the models share the structure of (Is]) and only differ in 
parameter values. Specifically, a true or nature run was 



created by integrating the Lorenz system 



xi = a{x2 


-Xi) 


X2 = rxi - 


X2 - X1X3 


is = ^1^2 


- bxs 



by direct insertion (see Fig. IT]). As an illustrative exam- 
ple, consider training with an analysis window of /i = 4. 
The first forecast, x^(to + 4), is a 4-timestep model fore- 
(8) cast started with the true initial condition x-^(to) = xq. 



with the standard parameter set: a = 10, b = 8/3, r = 
28. Models with the same a and 6, but with r- values 
varying from 25 to 31 in increments of 0.5 (except for 
r = 28) were the subjects for correction. For each of 
these 12 models, the correction algorithm was performed 
using the 4 different analysis windows h = 1,2,4, and 8 
timesteps, resulting in 48 distinct model-correction pairs 
in an exponential design. The training and testing of 
the corrected models is detailed in the following sections, 
and a picture showing one particularly positive outcome 
appears in Fig. |2j 




FIG. 2. (Color online) Trajectories of the truth (r = 28), 
a model with r-perturbation -2 (r = 26), and the corrected 
model using an analysis window of /i = 8 timesteps. All 
start from the same initial condition (circle), and represent 
2 time-units (200 timesteps) of integration (concluding with 
squares). Note that the corrected model trajectory is well 
aligned with the true trajectory for much longer than the un- 
corrected model trajectory. Even after it deviates noticeably, 
the corrected model trajectory changes flow regimes (switches 
lobes) with the true trajectory. In contrast, the uncorrected 
model trajectory deviates from the true trajectory almost im- 
mediately, and remains in the initial lobe. 



B. Training 

A 100 time-unit nature run was generated by in- 
tegrating system ([8| from the initial condition xq = 
[1.508870, -1.531271, 25.46091]^^ (see [18 ) for 10000 
timesteps of a>: = 0.01 time units each, using fourth order 
Runge-Kutta. For each of the 12 models and 4 analysis 
windows, a time-series of short forecasts was generated 



The second forecast, ^ {Iq 



is a 4-timestep model 



forecast started with the true state x^(to + 4), and so 
on, resulting in 10000/4 = 2500 total short forecasts. The 
state-independent correction {Ax.)/h was then computed 
as described in Fig. [l] and Sec. |IIB[ 

Next, the correction time-series Ax(t) was recomputed 
using a state-independent corrected version of the model, 
namely Eq. ([s]). Correction and true-state anomalies 
were calculated as in Eq. ([3|, and the cross covariance 
and covariance matrices were determined for each time t 
and averaged over the training sample to obtain C^xx^ 
and Gy.Ty.T as outlined in Sec. II C[ Finally, the static 
Leith operator L was computed and the training proce- 
dure was complete. 

Note that the training design imparts a statistical 
disadvantage upon the use of wider analysis windows. 
Specifically, doubling the analysis window halves the 
number of samples in the training period. The design was 
chosen, despite this prejudice, to more accurately reflect 
an operational implementation in which the training data 
is likely to be drawn from a fixed period of time. How- 
ever, to further support the validity of comparisons be- 
tween models corrected with different analysis windows, 
we note that letting the training period be lOOOO/i, en- 
suring that the number of samples is held constant at 
10000, yields results that are qualitatively indistinguish- 
able from those presented here. 



C. Testing 

A new nature run, 10000 time-units (one million 
timesteps) in length, was generated starting from the last 
true state in the training period. The purpose of begin- 
ning at the end of the training period was to obtain an 
independent truth with which to test the effectiveness of 
the corrected models. For each of the 48 corrected mod- 
els, 1000 randomly selected states from this new nature 
run were used as the initial state, and both the uncor- 
rected and corrected models were integrated for 20 time 
units. In addition, an exact model, with the same param- 
eters as the truth, was integrated from a random pertur- 
bation of that same initial condition, on the order of 10~^ 
in each state variable. The purpose of the exact model 
is to represent an upper bound for the effectiveness of 
the empirical model error correction, i.e. demonstrating 
a case where the impact of the r-perturbation was cor- 
rected exactly, and forecast error comes only from the 
initial condition discrepancy. 

Fig. [2] depicts the results of a single test with a par- 
ticularly positive outcome. Trajectories of the truth 
(solid black), uncorrected model (blue dot), and cor- 
rected model (red dash), all start from the same initial 
state and progress for two time-units. The corrected 
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FIG. 3. (Color online) Plots of average RMSE over 1000 tri- 
als for the uncorrected (thick solid blue) and corrected mod- 
els (all with r = 26) using analysis windows of 8, 4, 2 and 
1, black dot, red dash-dot, green dash, solid magenta respec- 
tively. The exact model (thin solid black) is the same as the 
true system, but started from an initial condition perturbed 
randomly on the order of 10~^ in each state variable. Time 
units are on the x-axis, where one time unit is 100 timesteps 
in the numerical integration. The technique's performance 
clearly degrades with widening analysis window. 



model trajectory stays close to the true trajectory for 
much longer than the uncorrected model trajectory, and 
even after deviating markedly, the corrected model tra- 
jectory still switches lobes with the true trajectory. 

Two metrics were used to measure forecast accuracy: 
root mean square error (RMSE), and anomaly correlation 
(AC). The RMSE is given at time t by 



RMSE(t) 



\ 



I E K(*) 



A^ 



a^fW]" 



(9) 



where x^ is the model state and x^ is the true state, 
and we are summing over the A^ = 3 state variables in 
the Lorenz system. Fig. [3] plots the average RMSE over 
1000 trials performed for correction of the r = 26 (r- 
perturbation of -2) model using the 4 different analysis 
windows. 

Anomaly correlation is a metric frequently used in 
weather and climate modeling to determine the length 
of time for which a model forecast is useful. The AC is 
given by 



,M' ^T' 



AC 



X- 



M' 



bllx^'lb 



(10) 



where x^ and x-^ are the anomalous model state and 
anomalous true state, respectively, at a particular time. 
AC is essentially the dot product of the anomalous model 
state with the anomalous true state, normalized such that 



AC = 1 for a perfect model. A forecast is typically con- 
sidered useful for as long as its AC remains above 0.6 [6j. 
As with RMSE, the AC scores for each model, corrected 
and uncorrected, were averaged over 1000 trials to pro- 
vide a good representation of model performance. See 
Fig. [4] for AC plots demonstrating the effects of changing 
analysis window length and parameter perturbation in 
the original model on the duration of useful forecasts. 
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FIG. 4. (Color online) (top) Plots of average AC over 1000 
trials for the uncorrected (thick solid blue) and corrected mod- 
els (all with r = 26) using analysis windows of /i = 8, 4, 2 
and 1, black dot, red dash-dot, green dash, solid magenta re- 
spectively. In agreement with Fig. [3] performance drops as 
the analysis window widens. Note, however, that even with 
/i = 4 timesteps the corrected model provides a useful forecast 
for twice as long as the uncorrected model, (bottom) Plots 
of average AC for models corrected with an analysis window 
of 1 timestep, and r = 25, 26 and 27, black dot, blue dash, 
and solid red respectively. The greater the magnitude of r- 
perturbation, the more difficult it is to correct successfully. 
Similar results are observed for r- values greater than 28. 



D. Results 



Duration of Useful Forecast 



Fig. [3] demonstrates that the tested empirical correc- 
tion technique succeeds in reducing error in the perfect 
model scenario. However, the importance of frequent 
observations of the truth during the training period is 
highlighted by the approach of the average corrected 
model RMSE towards the uncorrected model RMSE with 
widening of the analysis window. It is also noteworthy 
that the corrected model remains almost as good as the 
exact model for a full time unit when observations are 
made every timestep. In this case the RMSE remains 
below 2 more than 5 times as long the uncorrected model. 

Fig. [4] demonstrates that along with reduced error, em- 
pirical correction has the potential to provide forecasts 
that are useful for much longer. Training with an anal- 
ysis window of 1 timestep, the corrected model forecasts 
are useful for nearly 4 times longer than the uncorrected 
model forecasts. In light of the sensitivity of AC to anal- 
ysis window length, the bottom panel of Fig. |4] suggests 
that the accuracy of parameter values matters less for 
the effectiveness of the corrected model, as measured by 
error statistics, than does the frequency of observations 
in training. However, it should be noted that error statis- 
tics are not the whole story. The ability of the corrected 
models to reproduce the qualitative dynamical behav- 
ior of the true system accurately, like switching of lobes 
in the attractor, is only indirectly indicated by reduced 
forecast error. We address this issue in Sec. IIVDI 

A summary representation of the data in Fig. [4] is 
shown in Fig. [5J where the duration of forecast useful- 
ness is plotted vs. analysis window in the top panel, and 
vs. r-perturbation in the bottom panel. However, these 
plots still focus only on two cross sections of the 48 to- 
tal model-correction pairs that were tested in the per- 
fect model scenario, specifically the 4 corrected r = 26 
models in the top and the 12 corrections across all r- 
perturbations using an analysis window of 1 timestep in 
the bottom. To get a better sense of the relative impor- 
tance. Fig. [6] shows a surface plot of duration of useful 
forecasts for all combinations of analysis window and r- 
perturbation tested. Surprisingly, with an analysis win- 
dow of 2 timesteps, the corrected models with the largest 
r-perturbation (parameter error greater than 10%) pro- 
duce forecasts that are useful longer than those made by 
the uncorrected model with the smallest r-perturbation 
(parameter error less than 2%). For systems with reason- 
ably small model errors, this is an indication that empir- 
ical correction may improve forecasts more readily than 
parameter tuning. 



IV. TOY CLIMATE MODEL 

We next investigate the effectiveness of the correction 
procedure in a more realistic situation, where the forecast 
model is structurally different (in dynamics, dimension, 
parameterization, etc..) from the true system. Consider 
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FIG. 5. (Color online) Plots of the duration of useful fore- 
casts (first time AC=0.6) vs. analysis window (top) and r- 
perturbation (bottom), (top) The drop in duration of useful 
skill with analysis window suggests an asymptotic decrease 
in the value of correction towards no improvement from the 
uncorrected r = 26 model, (bottom) Using an analysis win- 
dow of 1 timestep, we plot the duration as a function of r- 
perturbation. The black diamond is the average duration of 
a useful forecast for the exact model (true system with per- 
turbed initial state), representing the limit of predictability 
imposed by initial condition uncertainty. Note: the corrected 
models with greatest r-perturbation outperform the uncor- 
rected models with least r-perturbation. 



a fluid-filled, vertically-oriented natural convection loop, 
or thermosyphon, with circular geometry, a schematic 
of which appears in Fig. [7[ The constant temperature 
imposed on the wall of the lower half of the loop, T/^, 
is greater than the constant temperature imposed on 
the wall of the upper half, Tc, resulting in a tempera- 
ture inversion. For large enough temperature differences 
AT^ = Th — Tc, convection dominates, and the flow un- 
dergoes chaotic reversals of direction referred to as Bow 
regime changes^ while remaining laminar [Ml |T5] . These 
dynamics produce forecasting difficulties similar to those 
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FIG. 6. (Color online) Surface plot showing duration of use- 
ful forecasts for all combinations of analysis window and r- 
perturbation tested, top surface corrected, bottom surface un- 
corrected. On the le ft face we can see the cross section repre- 
sented by Fig. 5(b) Using an analysis window of 2 timesteps. 



corrected models with largest r-perturbation provide forecasts 
that are useful longer than those of uncorrected models with 
smallest r-perturbation. 



encountered in weather and climate prediction, and thus 
the thermosyphon provides a useful platform on which 
to test potential improvements to forecasting methods. 




27r,0 



FIG. 7. Schematic of the natural convection loop, or ther- 
mosyphon. The top half of the fluid-filled loop is cooled to the 
temperature Tc^ while the bottom half is heated to Th > Tc, 
creating a temperature inversion. The three state variables 
in the low-dimensional model are (1) xi oc u, the mass flow 
rate (mean fluid velocity); (2) X2 oc ATs-g, the horizontal 
temperature difference; and (3) xs oc ATq-i2 — ATq1\^\ the 
difference between the vertical temperature profile and the 
value it takes during conduction. 



A. Experimental design 

The true system was represented by numerical simula- 
tion using the 2-D laminar Navier- Stokes equations along 
with the energy equation, and a finite- volume-based flow 
modeling software package (FLUENT 6.3) was used to 
perform the numerical integration, see [15] for details. 
Almost 90 days of fluid behavior was generated, in which 
O(IO^) flow reversals occurred. We also note that for the 
Rayleigh number used in this experiment, Ra = 1.5 x 10^, 
the thermosyphon has two unstable convective equilib- 
rium solutions corresponding to steady clockwise and 
steady counter-clockwise flow [151 [HI • 

The low-dimensional model used to make forecasts of 
the CFD simulation is the Ehrhard-Miiller (EM) system: 



Xi = a{x2 — Xi) 

±2 = j3xi- X2(l^ KH(\xi\)) - xix^ 

xs = xiX2-xs{l^KH{\xi\)) 



(11) 



where Xi is proportional to the mean fluid velocity, X2 is 
proportional to the horizontal temperature difference in 
the loop, and xs is proportional to the deviation of the 
vertical temperature difference from the value it takes 
during conduction. This system was derived from physi- 
cal principles to model a natural convection loop [TH [15] , 



and for this study the parameters a = 7, /3 = 33, K = 
0.07 were tuned empirically to best match the flow rever- 
sal behavior of the CFD simulated thermosyphon. 

The primary difference between EM and Lorenz sys- 
tems is H, a function that determines the velocity depen- 
dence of the heat transfer between the fluid and the wall. 
This characteristic of the flow is ignored by the Lorenz 
equations (i.e. K = 0). H varies as the third root of the 
magnitude of the mean fluid velocity for \xi\ > 1, and as 
a fourth degree polynomial in \xi\ for |xi| < 1 to remain 
different iable; the reader may see [15] for more detail. 



We note that when i^ = in the EM equations (11), 
they are identical to the Lorenz system ([8| with b = 1. 
Physically, the unitary geometric factor (i.e. b = 1) in 
EM results from the forced single circular convection cell 
in a thermosyphon, as opposed to the unconstrained flow 
producing multiple cells between two plates. 

Using a background forecast created with the EM 
model, and observations of the CFD mean fluid velocity 
u with Gaussian noise added to simulate error, 3DVar 
data assimilation was performed to generate an analysis, 
or best guess of the true state of the system in terms 
of the variables of the forecast model [15] . One segment 
of the analysis is used as the reference truth x-^(t) for 
training, and the remainder is used for testing. 



B. Training 

Approximately 3 days of 3DVar analysis, correspond- 
ing to 432 time-units in the EM forecast model, was used 
as the training period reference truth. As in the perfect 
model scenario, fourth-order Runge-Kutta was used with 
a timestep of k = 0.01 time-units. An analysis window 
of h = b timesteps, corresponding to about 30 seconds 
of simulated flow, was used to match the frequency of 
observation in the data assimilation scheme. Thus, 8640 
short forecasts were used to compute the model bias and 
Leith operator for the model, as outlined in Sec. |TI| 



C. Testing 

Three forecast models were compared by verification 
against both the 3DVar analysis and direct observation 
of the mass flow-rate u in the CFD simulated ther- 
mosyphon: (1) the uncorrected EM model with param- 
eters tuned to best represent observations of the CFD 
simulated mass flow rate; (2) the tuned model with cor- 
rection applied; and (3) an EM model whose parameters 
differ from the tuned model by 10%, with correction ap- 
plied. The purpose of the third test-model is to gauge 
the relative capabilities of empirical correction and pa- 
rameter tuning for error reduction and prolonging the 
usefulness of forecasts. 

A set of 1000 trial forecasts were performed, starting 
from randomly chosen points in the analysis after the 
training period. For each forecast, RMSE and AC time- 
series were computed with respect to the analysis and 
averaged over all trials. See the top panel of Fig. [S] for 
the resulting average AC plot (RMSE not shown). We 
also verify the model forecasts against observed scalar 
mass flow-rate u, for which time- series of relative error 
were averaged over the 1000 trials, pictured in the bottom 
panel of Fig. [8) 

Two details are important in the computation of the 
relative error. First, to compare model output to ob- 
servations, it is necessary to convert the model state- 
variables to "observation-space" variables. In other 
words, an observation operator determined by data as- 
similation was used to convert the model state-vector 
(xi,X2,X3)^^ to an observation-space value xi, which is 
the predicted mass flow-rate of the system. Second, the 
error is taken relative to the saturation point, which 
we deflne as the average absolute difference between the 
mass flow-rates of the system at randomly chosen points 
in time. Thus, an average relative error near 1 means that 
the forecast model is no better than a random guess. 

The results presented in Fig. [S] indicate that corrected 
models, tuned or not, produce smaller short-term fore- 
cast error on average than the uncorrected, optimally 
tuned model. Corrected model forecasts are thus typi- 
cally useful for longer. Though this is an important ben- 
eflt, average short-term error statistics may conceal con- 
siderable qualitative differences between model dynamics 
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FIG. 8. (Color online) Plots showing verification of the cor- 
rected model with respect to (top) 3DVar analysis and (bot- 
tom) direct observation of the CFD simulated mass fiow-rate 
u. (top) We see that empirical correction produces forecasts 
that are useful approximately 25% longer than forecasts of 
the uncorrected model. Also, correction of an untuned model 
(10% difference from optimal value in every parameter) pro- 
duces forecasts that are useful far longer than the uncorrected 
but optimally tuned model, (bottom) The corrected model 
shows reduced error when verified against direct observations 
of the mean fiuid velocity in the simulation as well. Again 
we see that correction of a model with sub-optimal parame- 
terization produces better results than tuning the parameters 
optimally but neglecting to apply empirical correction. 



and those underlying the true system. Stability of equi- 
librium solutions, and changes of flow regime character- 
ized by aperiodic switching between otherwise conflned 
regions of state-space are examples of qualitative char- 
acteristics for which it may be crucial that the model 
dynamics match those of the truth. In the next section 
we address the effect of empirical correction on the dy- 
namical matching capability of the EM forecast model. 



D. Matching qualitative dynamics 



To measure the accuracy of models with regard to 
matching the flow reversal behavior of the true system, 
forecasts were generated with both the corrected and 
uncorrected EM models for 5000 initial states through- 
out the attractor, from the testing portion of the 90- 
day 3DVar analysis, and the time of the first predicted 
flow reversal was recorded for each one. We investigate 
the difference between the predicted times and the ac- 
tual times (from the analysis) of first flow reversal, taken 
^modei — ^actuab SO that positivc valucs indicate late pre- 
dictions while negative values indicate early predictions. 
See Fig. [9] for plots of the results. 

Three unforeseen costs of empirical correction for this 
system, in terms of lost dynamical matching capability, 
can be summarized as follows: 

1) Stabilization of convective equilibrium solutions 

2) Elimination of flow reversal behavior for states in 
a neighborhood of either convective equilibrium 

3) Spurious dynamical asymmetry between lobes 
We address these costs in the sections that follow, exam- 



uncorrected model 



ining the first two related phenomena in Sec. |IVDT| an d 
then explaining the spurious asymmetry in Sec. |IVD2^ 
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1. Stabilization of equilibrium solutions 



Jacobian analysis of the EM equations (11) provides 



analytical confirmation of the instability of the convec- 
tive equilibria in the uncorrected model. Specifically, the 
Jacobian evaluated at each equilibrium has one negative 
real eigenvalue, whose eigenvectors are in the local direc- 
tion of (tangent to) the stable manifold of the equilib- 
rium, and a conjugate pair of complex eigenvalues with 
positive real part, whose 2-D eigenspace is locally tan- 
gent to the unstable manifold of the equilibrium. In fact, 
for both convective equilibria the positive real parts of 
these unstable eigenvalues are quite small, on the order 
of 10~^, indicating weakly repelling instability. In the fol- 
lowing we explain analytically the mechanism by which 
empirical correction overcomes this weak repulsion, pro- 
ducing a forecast model with attracting^ and thus stable. 



convective equilibria, see Fig. 10 



Empirical correction of the EM model effectively alters 
the right-hand side of the differential equations ( pT] ) by 
first adding a constant related to the bias term b, and 
then adding a term that depends linearly on the model 
state, i.e., something related to Lx^ Letting f^ be the 
vector-valued EM differential equation, and f^ be the 
corrected equation, we write 



where 



bo 



pM^ 



lim b 



pM 



bo + Lox' 



and 



Lo 



lim L 



(12) 



(13) 
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FIG. 9. (Color online) Difference (in minutes) between pre- 
dicted and actual time of first flow reversal, plotted by initial 
state, for the uncorrected (top) and corrected (bottom) mod- 
els. The difference was taken tmodei — tactual, so that positive 
values (towards red) indicate late predictions while negative 
values (towards blue) indicate early predictions. First, note 
that for initial states near the edge of the attractor, in either 
lobe, each model performs well in predicting the coming flow 
reversal. This makes sense because for these initial states the 
change of flow regimes is imminent, and thus easily predicted. 
Next, observe that when the uncorrected model errs, it almost 
always makes an early prediction (blue dots). Model trajec- 
tories close to the unstable convecting equilibrium solution at 
the center of each lobe oscillate with larger amplitude than 
in the true system. One might expect that empirical correc- 
tion would counteract the effect, and in fact that is apparent 
(bottom). However, rather than improving flow reversal pre- 
diction, the results here indicate extreme overcorrection. The 
dark red dots are not 40-minute late predictions; the limit of 
the color axis was chosen to better illustrate the spread in the 
remaining data. In actuality, the corrected model predicts 
that a flow reversal will never happen for those initial states. 
Finally, observe that the empirical correction has manufac- 
tured a lobe asymmetry, characterized by the much smaller 
region of dark red in the right lobe. 
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can be thought of as the computed bias term and Leith 
operator, respectively, for an infinitesimal timestep k. 
Because of the nonlinearity of the system, we cannot de- 
termine the exact relationship between the infinitesimal 
correctors bo and Lq, and the bias term b and Leith op- 
erator L, respectively, that we compute using a timestep 
of K, = 0.01 and analysis window of h = 5 timesteps. 
However, since we discretize f^ in numerical integra- 
tion, we can determine the bo and Lq that we actually 
apply. We effectively approximate the correction terms 



100b 



and 



lOOL (14) 



within the fourth order Runge-Kutta scheme. 

Now, armed with an analytical representation of the 
differential equations f^ , we note the following rela- 
tionship between the corrected model Jacobian Df^ , 
and that of the uncorrected EM model Df^: 



Dl 



M^ 



Di 



M 



Dl 



M 



lOOL 



(15) 



since the constant bias term bo disappears and Lo op- 
erates on a translation of the model state. We evaluate 
the Jacobian of the corrected model at each of the con- 
vective equilibria, and determine its eigenvalues. Indeed, 
the real part of the complex conjugate eigenvalues, for 
each convective equilibria, is negative for the corrected 
model f ^ , and thus the convective equilibria have be- 
come stable. 

The equilibria attract all states inside neighborhoods 
around them, which thereby separate state-space near 
the attractor into regions whose trajectories will either 
change fiow regime at least once, or not at all. See Fig. 
10 for an example of this effect. Furthermore, any trajec- 
tories that land in one of these neighborhoods after only 
one or two fiow reversals, which might occur within the 
expected 17- minute duration of useful forecasts for the 
corrected model reported in Fig. [S) will then approach 
steady convection. This behavior is in qualitative oppo- 
sition to that of the true system, and that of the original 
uncorrected EM model, for which steady convection in a 
single direction is an unstable equilibrium. 



In fact, empirical correction breaks this symmetry in two 
ways. 

First, recall that after bias correction alone we have 
changed the EM system by adding a constant vector 
to the right-hand side of the differential equation. Let- 
ting f^ represent the bias-corrected differential equa- 
tion, f ^ does not commute with A unless there is zero 
bias in xi and ^2, i.e. 



Af^*(x) ^ f^*(Ax) 



VboT^ 






^3 



(17) 



where 63 can be any constant, and we recall that bo = 
lim^^^o b can be thought of as the computed bias term 
for an infinitesimal timestep k,. Note that even if no bias 
in xi or X2 existed, the probability of statistically com- 
puting a bias term b that would preserve the symmetry 
of the EM system is zero. 

In the unlikely case that a bias term is computed that 
preserves symmetry, or such a bias term is forced, state- 
dependent correction will break it. Assuming that f^ 
does commute with A, and letting f^ be the fully cor- 
rected differential equation, then 



Af^^(x) 



f^^(Ax) 



?M+ 



ALo = LoA (18) 



In other words, f commutes with A if and only if Lo = 
lim^^o L commutes with A. This forces the computed L 
to be of the form 



L = 



hi 


^12 





^21 


^22 











^33] 



(19) 



where the £ij can be any constants. Of course the 
probability of computing such an L statistically is also 
zero. Therefore, both bias and state-dependent correc- 
tion break the symmetry of the EM model. 



3. Summary 



2. Broken symmetry 

The size discrepancy between left and right-lobe re- 
gions attracted to the convective equlibria of the cor- 
rected model revealed in Fig. |9] demonstrates that empir- 
ical correction breaks the symmetry of the EM system. 
As in the conventional Lorenz system (|8|, the EM sys- 
tem (11) is symmetric under the mapping (xi,X2,X3) \-^ 
(— xi, — X2,X3). Again letting f^ be the vector-valued 
EM differential equation, this symmetry implies that f ^ 
commutes with a certain matrix A, i.e. 



Ar^(x) = r^'(Ax), 



1 





0' 





-1 











1 



(16) 



Undesirable effects of empirical correction include the 
breaking of system symmetry, and altered stability of 
equilibrium solutions which results in large regions of 
state-space for which fiow-reversals never occur in cor- 
rected model forecasts. We emphasize that the qualita- 
tive behavior of the corrected model is substantively dif- 
ferent from that of the uncorrected model, which matches 
the behavior of the CFD simulated truth. This is true de- 
spite the fact that the corrected model shows improved 
average error statistics. However, it is possible to ad- 
just the correction procedure to mitigate this effect by 
directly incorporating dynamical knowledge of the true 
system, which is the subject of the next section. Note 
that in doing so, we sacrifice the general applicability of 
the technique. 
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1. Lobe dependence 

To generate lobe-dependent bias correction terms and 
state-dependent correction operators, two regions i^i,2 5 
corresponding to flow regimes of opposite direction, are 
defined by 



Li = {x : xi < 0}, 



{x : xi > 0} (20) 



noting that their union is the entire state-space. Physi- 
cally, I/i (L2) represents all states undergoing clockwise 
(counter-clockwise) convection. Next, the direct inser- 
tion procedure illustrated in Fig. [l] is modified to produce 
two subsequences of the analysis correction time series, 
Ax^^ and Ax^s, that correspond to two subsequences of 



the analysis time series, x^ 



and X? 

^2 



respectively, where 



FIG. 10. (Color online) The first 20 minutes of a forecast 
for an initial state attracted to the left-lobe convective equi- 
librium stabalized by the empirical correction (i.e., one of the 
dark red dots towards the center of the left lobe in the bot- 
tom of Fig. [9]). The uncorrected EM model trajectory (blue) 
behaves as it should, winding away from the equilibrium at 
the center of the lobe, whereas the trajectory of the corrected 
model (red) collapses toward the equilibrium solution, indi- 
cating the false stability produced by the correction. The 
source of the false stability is the over corrective nudging vis- 
ible as jumps (every 30 seconds, the length of the analysis 
window) on the red curve. The correction is attempting to 
lengthen flow regimes, a consequence of the blue regions near 
the equilibria in the top of Fig. [9] 



^Li 



{x^(t):x^(t)GLi} 



(21) 



and x|^ is defined similarly. These subsequences are sep- 
arated into mean and anomalous components, as in Eqs. 
(|3|, using means over each individual subsequence. Fi- 
nally, the separate correction terms and operators are 
computed by: 



h 



^Avii 



(22) 



for k = 1,2. To apply the lobe-dependent correction, 
at every timestep of numerical integration the current 
state is determined to be in either Li or I/2, and the 
appropriate bias correction term and state-dependent 
operator are applied to advance the model. 



E. Incorporating dynamical knowledge 



2. Dependence on distance from equilibrium 



To encode dynamical knowledge of the system in the 
empirical correction procedure, the state-space is parti- 
tioned into regions based on the qualitative behavior of 
the system, and then a separate bias term b and state- 
dependent correction operator L is computed for each 
region. For example, in the context of weather forecast- 
ing, the state-space of the atmosphere could be divided 
by stage in the El Nino oscillation, or by day and night, or 
local season for regional models. Fig. |9] suggests two ways 
to partition the state-space in the present context: (1) 
by flow regime direction (lobe); or (2) by distance from 
the nearest convective equilibrium solution. In each case, 
the state-space is decomposed into two regions, left/right 
lobe, or near/far from equilibrium, respectively. In ad- 
dition to testing the correction procedure using each of 
these strategies individually, a procedure applying them 
simultaneously, which results in a partition of the state- 
space into four regions, is also tested. 



A procedure analogous to that for lobe-dependent cor- 
rection is applied here. Two regions £^1,2, corresponding 
to near and far from equilibrium, respectively, are defined 

by 



El = {x : min(||x — eqi| 
E2 = {x:x^Ei} 



■eq2||2. 



< Cy}, 



(23) 



where eqi and eq2 are the convective equilibria (esti- 
mated from the parameters of uncorrected model, [H]), 
and the critical value Cy is a parameter for the procedure. 
We are effectively approximating the neighborhoods at- 
tracted to the convective equilibria by spheres of radius 
Cy. For all results shown in the paper, the critical value 
Cy = 8.5 was used, though error statistics and dynamical 
matching capability were virtually unaffected by chang- 
ing this parameter by 25% in either direction. This range 
of critical values was tested as estimates to the average 
radius of the dark red region in the left lobe of Fig. [9] 
(bottom). Continuing with the correction scheme, direct 
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insertion is modified as in the lobe-dependent correction, 
and the region-specific bias terms and state-dependent 
operators are calculated as in Eqs. (22), substituting Ej^ 



for L/c. Application of the correction to a forecast model 
also proceeds in the same fashion. 



3. Simultaneous lobe and equilibrium dependence 

Defining the lobe regions Li^2 as in the lobe-dependent 
section, and the equilibrium regions £^1^2 as in the previ- 
ous section, we define the four regions for simultaneous 
lobe and equilibrium-dependent correction by 



4 = LiDEi, 


i?2 


= LinE2 


,3 = L2nEi, 


Ri 


= L2 n £^2 



(24) 



SO that we modify direct insertion to produce four subse- 
quences of analysis increments, each paired with the ap- 
propriate subsequence of the analysis time series. Note 
that the critical value Cy used in defining the regions £^1^2 
does not depend on the lobe in this scheme. Allowing a 
different Cy for each lobe is a possible modification that 
was not tested. We compute the bias terms and Leith op- 
erators as in Eqs. (22 ), substituting Rk for L^, where now 
k = 1,2,3,4. Again, application of the correction pro- 
ceeds as in the individual dynamically informed schemes, 
where the current state is determined to be in one of the 
four defined regions, and the appropriate bias term and 
Leith operator are used to advance the model. 



perfect model scenario. Verification against both analy- 
sis and observations allows more confidence in the suc- 
cess of the procedure in aligning forecasts made by the 
low-dimensional model with the CFD simulated ther- 
mosyphon. 

Additionally, the results constitute evidence that the 
correction procedure is more effective than fine-tuning of 
parameters for improving error statistics in this more re- 
alistic setting, as well as in the perfect model scenario. 
This is not to say that parameters should not be tuned, 
but rather that empirical correction is a cheaper and 
more effective avenue for reducing forecast uncertainty 
in the present context. 

Empirical correction of the EM forecast model comes 
at the cost of decoupling qualitative dynamics from those 
of the true system, however. As shown in Fig. |9J a large 
region is created in which initial states have trajectories 
under the forecast model that behave completely differ- 
ent from how they would in the true system. An example 
of such a trajectory appears in Fig. 10 While the gen- 



eral empirical correction procedure reduces forecast error 
for many initial states, it produces a forecast model that 
is entirely useless for others. By introducing dynamical 
knowledge into the correction procedure, the results of 
which appear in the next section, we greatly reduce the 
number of such initial states, and also further reduce av- 
erage forecast error. 



2. Dynamically informed correction 



F. Results 

First, we discuss the results of applying the dynam- 
ically uninformed, and entirely general correction pro- 
cedure as described in Sec. [n[ to couple the three- 
dimensional EM forecast model to the high-dimensional 
CFD simulated thermosyphon. Then we report the re- 
sults of modifying the correction procedure to directly 
incorporate dynamical knowledge of the system, as de- 
tailed in Sec. ITVEl 



1 . General procedure 

The ability of the correction procedure to overcome 
parameter inaccuracies, as suggested by the results in 
the perfect model scenario, inspired the comparison of 
three forecast models. 

(1) tuned, uncorrected 

(2) untuned, corrected 

(3) tuned and corrected 

Fig. [8] shows that empirical correction produces forecasts 
that remain useful longer, and demonstrate reduced er- 
ror in this mock-operational setting, as well as in the 



In an effort to produce forecast models that more accu- 
rately refiect the qualitative dynamics of the true system, 
two types of system-specific dynamical knowledge were 
directly incorporated into the empirical correction pro- 
cedure: (1) flow regime direction of the current system- 
state, and (2) distance between the state and the nearest 
convective equilibrium solution. These dynamical cues 
were incorporated individually, and also simultaneously, 
resulting in three dynamically informed, empirically cor- 
rected forecast models, as outlined in Sec. |IVE[ 

Fig. 11 shows average AC and mass flow-rate relative 
error over 5000 trials for the three dynamically informed 
and corrected models, as compared to the original bi- 
ased model and the dynamically uninformed, corrected 
model. As described in the caption to the flgure, encod- 
ing the current flow regime into the correction procedure 
results in a forecast model that is useful for almost twice 
as long as the original biased model, and doubles the im- 
provement that was gained by dynamically uninformed 
emprical correction. Encoding the distance to the nearest 
equilibrium solution, on the other hand, does not greatly 
prolong usefulness, and in fact reduces it slightly when 
applied simultaneously with lobe-dependent correction. 

In Fig.[T2J we see how the three forecast models result- 
ing from dynamically informed empirical correction com- 
pare with the dynamically uninformed, corrected model, 
with regard to matching the qualitative behavior of the 
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FIG. 11. (Color online) Average AC (top) and mass flow- 
rate relative error (bottom) over 5000 trials for models that 
were: uncorrected (thick solid blue), corrected (green dash- 
dot), equilibrium-dependent (ED) corrected (cyan dash), lobe 
and equilibrium-dependent (LD-ED) corrected (black dot), 
and lobe-dependent (LD) corrected (solid red). Note that the 
models are listed in order of increasing performance, as mea- 
sured by these error statistics. Lobe-dependent correction, 
i.e., correcting with bias terms and state-dependent opera- 
tors that depend on the current flow regime direction, pro- 
duces forecasts that are useful for almost twice as long as 
those made by the uncorrected model, on average. When 
compared to the dynamically uninformed correction scheme 
(green dash-dot), lobe-dependent correction more than dou- 
bles the improvement over the uncorrected model. We note 
that equilibrium-dependent correction does not produce much 
improvement in forecast error over the dynamically unin- 
formed procedure, and hinders the performance when applied 
simultaneously with lobe-dependent correction. 



true system. We see a trend of improvement, charac- 
terized by smaller regions of states whose dynamics are 
different in the models than in the CFD simulated ther- 
mosyphon (dark red regions), as we apply first lobe- 
dependent, then equilibrium-dependent, and finally si- 
multaneously lobe and equilibrium-dependent correction. 



dure primarily enhances forecast statistics, while encod- 

TABLE L Median absolute differences between predicted 
and actual times of flrst flow reversal (row 1), along with 
percentage of trials for which the flrst flow reversal was pre- 
dicted within 1 and 2 minutes, (rows 2 and 3, respectively), 
for the uncorrected model (M), dynamically uninformed cor- 
rected model (CM), lobe-dependent corrected model (LD), 
equilibrium-dependent corrected model (ED), and simultane- 
ously lobe and equilibrium-dependent corrected model (LD- 
ED). 





M 


CM 


LD 


ED 


LD-ED 


median 


4 


10 


1.5 


10.5 


1.5 


% within Imin. 


38.26 


38.54 


49.12 


38.68 


48.64 


% within 2min. 


47.04 


43.40 


54.38 


43.58 


57.08 



ing distance from equilibrium primarily enhances dynam- 
ical matching. Simultaneous inclusion of the two types 
of dynamical cues results in the best dynamical match- 
ing, with only a slight cost in average forecast accuracy. 
For further evidence of this summary conclusion, consider 
Table [D 

In the first row, the median absolute differences 
between predicted and actual times of first flow re- 
versal over the 5000 trials are listed for the uncor- 
rected model (M), dynamically uninformed corrected 
model (CM), lobe-dependent corrected model (LD), 
equilibrium-dependent corrected model (ED) and simul- 
taneously lobe and equilibrium-dependent model (LD- 
ED). The means are highly skewed for the CM, LD and 
ED models, due primarily to the large basins of attrac- 



tion for the left-lobe equilibrium solution, see Fig. [12] top 
left, right and bottom left, respectively, and thus we show 
only the medians. The medians show that both the LD 
and the LD-ED corrected models predict the first flow re- 
versal more accurately than the uncorrected model. The 
bottom two rows show the percentage of the 5000 trials 
for which the first flow reversal was predicted within 1 
and 2 minutes, respectively, for each of the models. Again 
the LD and LD-ED corrected models show the best per- 
formance. Note that the LD-ED model boosts both the 
1 and 2-minute success rates by approximately 10% over 
the uncorrected model. 

We note, however, that even the LD-ED model ex- 
hibits a spuriously stable convective equilibria in each 
lobe. Fig. 12 Incorporating dynamical knowledge of the 
true system in the correction procedure, at least through 
the partitioning of state-space as we have done here, is 
not enough to avoid the stabilizing effect of empirical 
correction on the equilibria of the EM model explained 
in Sec. 



IVDl 



Results shown in Figures 11 and 12 suggest that en- 
coding ffow regime direction into the correction proce- 



It is plausible that models of other sys- 
tems with weakly repelling (attracting) equilibria might 
be subject to similar stabilizing (destabilizing) effects un- 
der empirical correction. Such effects might be mitigated 
by the strategy employed in the present work, but prob- 
ably not avoided. 
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FIG. 12. (Color online) Difference (in minutes) between predicted and actual time of first flow reversal, plotted by initial 
state, for the (A) corrected, (B) lobe-dependent corrected, (C) equilibrium-dependent, and (D) lobe and equilibrium-dependent 
models. As in Fig. ^ the difference was taken tmodei — tactual, so that positive values (towards red) indicate late predictions 
while negative values (towards blue) indicate early predictions. Inset partial histograms show the number of forecasts (out 
of 5000) predicting the first flow reversal within 3 minutes of the truth (green) and predicting that a flow reversal will never 
happen (red); the axes have the same scale for comparison, and the bar colors correspond to the dot colors. (B) Applying 
lobe-dependent correction (as opposed to (A) dynamically uninformed correction) increases the number of initial states for 
which the first flow reversal is predicted accurately, visible as a larger number of green dots (bigger green bar). Also, although 
the region of dark red (initial states attracted to convective equilibrium) in the left lobe has decreased in size, the one in the 
right lobe has inflated, again as compared to the top left, (bottom left) Equilibrium-dependent correction shrinks the left-lobe 
red region without inflating the one in the right lobe. However, it maintains the region of initial states for which flow reversal 
predictions are slightly early (light blue), which was reduced by the lobe-dependent correction, (bottom right) Applying the 
simultaneously lobe and equilibrium dependent correction, in which there are four different correction regions, the forecast 
model demonstrates the smallest region of initial states whose qualitative dynamics are different from the CFD simulated 
thermosyphon. 
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CONCLUDING DISCUSSION 



It is apparent from the results of this work that the 
empirical correction technique tested here is successful 
in reducing forecast error in the low-dimensional set- 
ting. In both the perfect model scenario and the mock- 
operational experiment, improved error statistics and 
prolonged usefulness of forecasts were demonstrated by 
the corrected models. Furthermore, the empirical cor- 



rection procedure was shown to provide greater improve- 
ment in average forecast accuracy than fine-tuning of pa- 
rameters. However, this positive result comes at the cost 
of altering some important dynamical characteristics of 
the model. Details and implications of these results are 
discussed in the next two sections. In the third sec- 
tion, we make an observation about the relative impact 
of state-independent and state-dependent correction, in 
the present context, to conclude the paper. 
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A. Empirical correction vs. parameter tuning 

In the perfect model scenario, empirical correction of 
models with the greatest parameter error provided fore- 
casts with smaller short-term average error, and that 
were useful for longer, than those made by the uncor- 
rected models with least parameter error. Similarly, em- 
pirically correcting the EM model with 10% error in every 
parameter (measured from tuned values) produced a fore- 
cast model that outperforms the tuned, but uncorrected 
model in predicting the CFD simulated thermosyphon. 
In each experiment, superior performance is evident in 
terms of anomaly correlation and average forecast error, 
measured against 3DVar analysis, and is verified against 
observed mass fiow-rate of the CFD simulation in the toy 
climate experiment. 

We do not present these results as evidence that em- 
pirical correction should replace parameter tuning, nor 
even that the former is better than the latter in any well- 
defined sense. However, the results do suggest that em- 
pirical correction could be a viable complement to the 
tuning of model parameters. Particularly as degrees of 
freedom become large, for example N ^ 10^^ in some cur- 
rently operational numerical weather models, the compu- 
tational cost of parameter tuning is very large in compar- 
ison to that of empirical correction, when appropriately 
modified for such models (see O US]). An example of 
such a modification is to compute the first m <^ N prin- 
cipal components (PCs) of L by singular value decom- 
position (SVD), where m is determined so that a certain 
percentage of the state covariance is explained by these 
first m PCs. A combined tuning/correction approach 
could reduce the number of model integrations necessary 
in the parameter tuning process without sacrificing model 
accuracy. However, the dynamical ramifications of this 
strategy must be considered, as we explain in the next 
section. 



B. Error statistics vs. system dynamics 

The reduction in average forecast error provided by 
empirical correction belies fundamental dynamical dis- 
turbances born out of the correction procedure. Stabi- 
lization of equilibrium solutions, which results in large 
regions of initial states for which the corrected model 
forecasts behavior in opposition to that of the true sys- 
tem, follows from the dynamical modification imposed by 
empirical correction. In addition, the symmetry of the 
model system is broken by empirical correction. Though 
these costs can be mitigated somewhat by hard-wiring 
system-specific dynamical cues into the correction pro- 
cedure, they can not be eradicated without more funda- 
mental alterations of the technique, e.g., forcing the bias 
term b and Leith operator L to preserve system sym- 
metry. In operational practice, empirical correction is 
known to introduce imbalances, e.g. violating geostro- 
phy, necessitating some mechanism for smoothing the 



fiow into a physically viable region of state space. 

In fact, it may be impossible to avoid all dynamical in- 
accuracies resulting from empirical correction, and even 
if theoretically possible, it would likely be impractical to 
do so in any operational setting. In considering the ap- 
plication of the technique in operational settings, then, 
it must be determined if the effects of misrepresented dy- 
namics can be reduced to a tolerable level on a case-by- 
case basis. In the ideal situation, regions of state-space 
that would be dynamically misrepresented under empir- 
ical correction could be reduced, by minor modifications 
to the correction procedure, to encompass only unrealis- 
tic or unlikely physical states. In any case, the technique 
presented in this study should not be applied without 
such considerations. 



C. Bias correction vs. Leith operator 

State-independent error correction by itself produces 
almost no improvement in any of the forecast models 
in this study (not shown). This is in contrast to what 
has been observed in operational weather and climate 
model studies, where state-independent bias correction 
typically outperforms state-dependent correction in re- 
duction of forecast errors [20]. The inaccuracies of ad- 
hoc forcings included in such models to compensate for 
external and/or irresolvable phenomena (e.g. solar and 
cloud forcings, respectively) are likely responsible for a 
large component of the bias. In light of the lack, or at 
least minimal nature of such external and sub-gridscale 
infiuences in the toy models considered here, the ineffec- 
tiveness of bias correction is logically consistent with this 
explanation. 

The state-dependent Leith operator is entirely respon- 
sible for the success of the corrected models in this 
study. In the perfect model scenario this makes sense 
because the difference between the forecast models and 
the "truth" model are inherently multiplicative, i.e. the 
parameters are coefficients weighting the interaction be- 
tween state- variable values and thus resulting errors must 
depend on state. For the EM model of the CFD system, 
it seems that errors resulting from the low dimensional- 
ity of the forecast model may also be multiplicative in 
nature. If this is the case, state-dependent correction 
may reduce error patterns in operational models that re- 
sult from reduced dimensionality, e.g. coarse resolution. 
The correction will not likely compensate for processes 
that are irresolvable due to coarse resolution, but rather 
may reduce the propagation of error resulting from the 
omission of such phenomena. This hypothesis is consis- 
tent with demonstrated improvement of local behavior 
in state-dependent corrected atmospheric models with 
N ^ 10^ degrees of freedom |9.. 

In previous studies of state-dependent correction in 
models that are much more realistic than those consid- 
ered here, resulting error reduction has been minuscule in 
comparison to what is achieved by bias correction. How- 
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ever, this is not cause to reject the usefulness of parame- 
terizing state-dependent error. Though globally averaged 
error reduction may not be significant, improvement in 
the local behavior of models can have a large impact on 
forecast uncertainty, particularly in an ensemble strategy 
where state-dependent correction can increase the spread 
in previously unsampled state-space directions. 
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